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Abstract In time-distance hclioscismology, wave travel times arc measured from 
the cross-correlation between Dopplcr velocities recorded at any two locations 
on the solar surface. However, one of the main uncertainties associated with 
such measurements is how to interpret observations made in regions of strong 
magnetic field. Isolating the effects of the magnetic field from thermal or sound- 
speed perturbations has proved to be quite complex and has yet to yield reliable 
results when extracting travel times from the cross-correlation function. One 
possible way to decouple these effects is by using a 3D sunspot model based 
on observed surface magnetic-field profiles, with a surrounding stratified, quiet- 
Sun atmosphere to model the magneto-acoustic ray propagation, and analyze 
the resulting ray travel-time perturbations that will directly account for wave- 
speed variations produced by the magnetic field. These artificial travel-time 
perturbation profiles provide us with several related but distinct observations: 
i) that strong surface magnetic fields have a dual effect on hclioseismic rays - 
increasing their skip distance while at the same time speeding them up con- 
siderably compared to their quiet- Sun counterparts, ii) there is a clear and 
significant frequency dependence of both skip-distance and travel-time pertur- 
bations across the simulated sunspot radius, in) the negative sign and magnitude 
of these perturbations appears to be directly related to the sunspot magnetic- 
field strength and inclination, iv) by "switching off" the magnetic field inside 
the sunspot, we are able to completely isolate the thermal component of the 
travel-time perturbations observed, which is seen to be both opposite in sign 
and much smaller in magnitude than those measured when the magnetic field 
is present. These results tend to suggest that purely thermal perturbations are 
unlikely to be the main effect seen in travel times through sunspots and that 
strong, near-surface magnetic fields may be directly and significantly altering 
the magnitude and lateral extent of sound-speed inversions of sunspots made by 
time-distance helioseismology. 
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1. Introduction 

Time-distance helioseismology is a powerful diagnostic tool used in local hclio- 
seismology to probe the subsurface structure and dynamics of the solar interior, 
in particular in and around solar active regions. To date however, results obtained 
by time-distance helioseismology have not directly accounted for the effects of the 
magnetic field on the wave-speed in travel-time perturbation maps, forward mod- 
elling or inversions, but have indirectly included magnetic effects only through 
their influence on the acoustic properties of the medium (e.g. the sound speed). 
Standard forward-modelling is based on a number of assumptions including, but 
not limited to, Fermat's Principle and the ray approximation (e.g. Kosovichev, 
Duvall, and Scherrcr, 2000; Zhao, Kosovichev, and Duvall, 2001; Hughes, Ra- 
jaguru, and Thompson, 2005), the Fresnel-Zonc approximation (e.g. Jensen et 
at, 2001; Couvidat et ah, 2004) and the Born approximation (e.g. Couvidat, 
Birch, and Kosovichev, 2006). These models do not include any provision for 
surface effects. In fact, no standard local-helioseismic method includes provisions 
for contributions from near-surface magnetic fields. 

Recent work in sunspot seismology has pointed to the significant influence of 
near-surface magnetic fields and possible contamination due to their effects in hc- 
lioseismic inversions for sound speed beneath sunspots (Couvidat and Rajaguru, 
2007). Prior to this, a number of other very important results have highlighted 
the complications of interpreting helioseismic observations (in particular, the 
interaction of p modes) in the near-surface regions of sunspots (see e.g. Fan, 
Braun, and Chou, 1995; Cally, Crouch, and Braun, 2003; Lindsey and Braun, 
2005; Schunker et ai, 2005; Schunker and Cally, 2006; Braun and Birch, 2006). 

The key issues are i) how to successfully model the effects of wave-speed 
inhomogcncitics thought to be produced by the magnetic field in solar active 
regions, ii) how to isolate such effects from those thought to be associated with 
temperature, flow perturbations, and other observational constraints and effects, 
and finally Hi) how will inferences made about subsurface structure change as 
a result of incorporating these effects into the modelling process? Efforts to 
address these issues both observationally and computationally have been largely 
unsuccessful, mainly because of a general lack of understanding of the process 
involved. But there is some light at the end of the tunnel, as there are cur- 
rently under development a number of robust magnetohydrodynamical (MHD) 
simulations modelling helioseismic data and wave propagation that may aid our 
understanding considerably in the near future (e.g. Cameron, Gizon, and Duvall, 
2008; Hanasoge and Duvall, 2007). In this work, we shall attempt to address 
some of these outstanding issues by using helioseismic ray theory to perform 
forward modelling of helioseismic rays in a simulated sunspot atmosphere with 
the aim of modelling the magneto-acoustic ray propagation and analysing the 
resulting artificial ray travel-time perturbations that will directly account for 
wave-speed variations produced by the magnetic field. We shall also address 
the problem of trying to isolate and analyze the thermal contributions to the 
observed travel-time perturbations using our simulations. 
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2. The Sunspot Model 

The axisymmetric sunspot model chosen for our analysis consists of a non- 
potential, untwisted, magnetohydrostatic sunspot model constrained to fit ob- 
served surface magnetic field profiles. The surface field is therefore quite realistic, 
which is important because there is evidence (Schunker and Cally, 2006) that 
magnetic effects on helioseismology are dominated by the top Mm. 
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Figure 1. Plots of the radial (B r , left), vertical components of the observed magnetic field 
(B z , right) and magnetic field inclination from the vertical (8°) as derived from IVM surface 
magnetic field profiles of Active Region (AR) 9026 on 5 June 2000, shown as a function of 
sunspot radius (r, Mm). Solid lines indicate constrained polynomial fits. Values of B are shown 
in Gauss (G). 

The sunspot also needs to be surrounded by an unperturbed, stratified atmo- 
sphere. The background model employed consists of a Global Oscillation Network 
Group (GONG) Model S atmosphere (Christensen-Dalsgaard et al, 1996) (ob- 
tained from the {L5BI.D.15C.PRES.960126.AARHUS} Model S package). The 
preferred surface field configuration of the flux tube was derived from constrained 
polynomial fits to the observed scatter plots of the radial (B r ) and vertical (B z ) 
surface magnetic field profiles (see Figure 1) of AR 9026 on 5 June 2000 - a fairly 
symmetrical sunspot near disk-centre, ideal for hclioseismic analysis - obtained 
from IVM (Imaging Vector Magnetograph) vector magnetograms (see Mickey 
et al. (1996) for more details regarding the observations). We note that in the 
B z profile of AR 9026, the vertical-field strength tends to decrease to around 
2 kG as it approaches r = 0. We find this highly improbable for a sunspot, 
so we extrapolate to a peak field of 3 kG for our model at r = 0. (A separate 
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analysis was conducted for the model with the (unrealistic) peak field of 2 kG. As 
expected, the only difference we observed was the magnitude of the perturbations 
produced being slightly smaller than the ones we report in Section 4. All other 
results appeared to be identical) . The fits of B r and B z are then used to derive 
an analytical form for the potential function, 

^(r,z)=^o(^fr) (!) 



n(z) 

where Vo is the derived surface field at the surface (z = Z ), the radius of the 
sunspot at the surface (r = R ) is fixed at Ro = 16 Mm. Instead of a current 
sheet along the boundary, we prescribe an analytical form for the outermost field 
line, 

n(z) = f° r z R ?\i\ i hi? ™> ( 2 ) 

where the field strength drops to zero and R m and c are free parametres. We 
ensure that all calculations (e.g. change in pressure, density, etc.) made across 
the boundary layer/transition region between the sunspot atmosphere and the 
external environment are both consistent and continuous along r^. 

The next step essentially involves solving the standard equations of magnc- 
tohydrostatics (MHS), using the Model S atmosphere and its variables as the 
quiet-Sun environment. The magnetic pressure and tension resulting from the 
Lorentz force, 

f L = J x B, (3) 

are confined within the simulated sunspot atmosphere, where p is the magnetic 
permeability and J = jj(V x B). The gas pressure p(r,z) is calculated using 
horizontal force balance, 

Pi{r,z) =p e {z) + Ap(r, z) (4) 

where pi(r,z) and p e (z) denote internal and external (i.e. Model S) pressure 
respectively and the change in pressure is therefore 

Ap(r,z) = f f Lr dr (5) 

which drops to zero as we approach r b . Once the pressure inside the sunspot and 
along the boundary are known, the density p(r, z) , can similarly be calculated 
using vertical force balance, 

Pi(r,z) = p e (z) + Ap(r,z) (6) 

where the change in density is given by 



Ap(r, z )= l - 



' dAp(r,z) 

JL Z - 



dz 



(7) 
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Figure 2. The magnetic field configuration for the sunspot model. The field lines plotted indi- 
cate equidistant magnetic-flux values. Internal and external (Model S) variables are indicated 
for reference. r\> represents the radius of the outermost field line, which varies with depth (z) 
along the sunspot radius. 



This is essentially all that is required to then compute the modified sound 
speed or thermal profile of the sunspot atmosphere, 



c 2 s (r,z) = c 2 s (z) + T 1 (z) 



Pi(r,z) _ p e (z) 
Pi(r,z) p e (z) 



(8) 



while for the sake of simplicity, assuming the ratio of specific heat (Ti) that 
appears in the sound speed is the same function of height as it is in the external 
atmosphere. Finally, all that is left is to calculate the Alfven speed, 



a 2 (r, z) 



1 



[B 2 + B 2 ]. 



(9) 



Some of the important internal properties of the resulting sunspot model (e.g. 
pressure, density, sound and Alfven speeds) are shown in Figure 3. The external 
(Model S) profiles for each variable are also shown for reference. The near-surface 
thermal structure of the sunspot and the (a = c s ) equipartition depth is also 
shown for reference in Figure 4. We can clearly see the modified sound-speed 
structure (c 2 ) as a result of the magnetic field in this image. It is interesting 
to note that in our (simple) model the region of decreased sound-speed does 
not appear to extend as deep as 3D time-distance inversions of the real Sun 
have suggested. Estimates for the lateral extent of the decreased sound-speed 
region using tomographic imaging of the sub-surface layers of sunspots have 
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ranged from depths of approximately z = —2.4 to z = —3.5 Mm using the Born 
and ray approximations respectively (Couvidat, Birch, and Kosovichev, 2006). 
Nevertheless, the sunspot model exhibits the broad features expected of a real 
sunspot, and presents a useful test case. 
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Figure 3. Internal pressure (p), density (p), sound (c s ), and Alfven (a) speed profiles of the 
sunspot model with an external GONG Model S atmosphere. Left-hand coloumn profiles are 
calculated along the surface of the sunspot (z = 0), while right-hand coloumn profiles are 
calculated along the axis of the sunspot (r = 0). Internal profiles are indicated by solid lines 
in all plots. The thick solid line in the bottom two panels indicate Alfven speeds. The dashed 
lines represent GONG Model S values in all plots. 



3. Ray Path Calculations 

The ray paths are calculated in Cartesian geometry, in the realm of frequency 
dependent ray paths described by Barnes and Cally (2001), with the complete 
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Figure 4. The thermal profile (c^) in the top 1 Mm of the sunspot. Lighter coloured contours 
(i.e. cyan/green) indicate regions of decreased sound speed (cooler regions) under the sunspot 
surface, while darker (hotter) regions (i.e. orange/red) are indicative of areas of enhanced sound 
speed. The dashed line marks the position of the a = c s layer. Field lines are over-plotted 



form of the three-dimensional dispersion relation: 

V = u 2 u 2 c a 2 y k 2 h + {u 2 - a 2 k\) x [co 4 - (a 2 + c 2 )u 2 k 2 

+a 2 c 2 k 2 kf + c 2 N 2 k 2 h - (oj 2 - a 2 z k 2 )oj 2 c ] = 0, (10) 

where kh and fcii are the horizontal and parallel components of the wave-vector 
k and 

N 2 = j-- 9 - 2 (11) 
Hp & 

is the squared Brunt- Vaisala frequency, with g being the gravitational accelera- 
tion, Hp(z) the density scale height, and H' p = dH p /dz and ui 2 is the square of 
the acoustic-cutoff frequency. For completeness, we calculate the raypaths using 
two forms of u> c . The most commonly used form 



.2 



(1-2F1), (12) 



exhibits an extended sharp spike around z = —100 km (see Figure 5). This form 
of ui c is often used by helioseismologists. However, as Cally (2007) points out, 
this sharp spike in the cutoff frequency is inconsistent with the WKB assumption 
of slowly varying coefficients on which T> is based. A much smoother isothermal 
form, 

w c . = c/2H, (13) 

is consistent with the derivation of T>, and does not suffer from the spike (see 
Figure 5). Unless otherwise stated, all results shown here utilize u Ci . (Simula- 
tions using the form of lo c in Equation (12) were also conducted, the results 
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Figure 5. Plots of the various forms of the acoustic cutoff (lu c ) and Brunt-Vaisiila (iV) 
frequencies. The later is indicated by a blue, solid line inside the sunspot atmosphere and 
dashed blue line indicating Model S values. The solid black line indicates the acoustic cutoff 
frequency u) c for the sunspot atmosphere, while the dashed black line indicates Model S values. 
The isothermal form, ui Ci is indicated by the solid red line for the sunspot atmosphere, dashed 
red line indicates Model S values. 



being very similar to those reported in Section 4, expect for a certain amount 
of unsmoothness being present in the travel-time perturbation profiles (mainly 
affecting shallow rays which are more sensitive to the reflecting boundary near 
the surface) as a result of using the more rigid form of lo c ). Naturally, the 
magnetic field slightly modifies both uj c and lu c . , the results of which can be 
seen in Figure 5. 

Following Weinberg (1962), the construction of k is completed by specifying 
the governing equations of the ray paths 



dx _ dV 
d7 ~ dk 



(14) 



dk 


&D 


dr ~ 


dx 


dt 


dV 


d7 ~ 


du 


duj 


dV 


d7 = 





(15) 
(16) 
(17) 

where r parameterizes the progress of a disturbance along the ray path. For a 
time-independent medium, for which dT>/dt = and lu is constant, the phase 
function S(x) evolves according to 

dS 1 , dx . . 

— =k-— -u. (18) 
dt dt y ' 
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Hence, 

S(x) = J k-dx-wi, (19) 

where the first term (integral) represents the contribution to the phase due to 
motion along the ray path, and the second term represents the Eulerian part. 
Since we are only going to be concerned about the change in phase due to motion 
along the ray path, we can essentially ignore the Eulerian part for the rest of 
our analysis. 

We iteratively find the initial wave- vector (ki n ; t ) by using an initial guess 
which comes from solving T> = for the wavenumber, assuming the wavevector 
is in the directions a, (3 - where a and j3 are angles from the vertical and the 
x-z plane respectively of the initial shot. Initially, we initiated the rays from the 
top of the ray path, adjusting the initial shooting angle (a) to obtain the desired 
range of skip distances. However, given the very sensitive nature of the near- 
surface region of the sunspot atmosphere, we used a much finer computational 
grid in the top 1.5 Mm. As a result, we encountered many instances of rays 
initiated inside evanescent regions (which should obviously be avoided) and also 
obtaining very shallow rays with little or no helioscismic value. So in order 
to reduce computation time and also have greater flexibility in choosing the 
desired range of ray skip distances, we initialized the rays from the minima of 
their trajectories (essentially the lower turning point of the ray, z bot ). Hence, 
the value of a was fixed at a = 90°, allowing us to adjust the initial shooting 
depth Zbot to obtain the desired range of skip distances. 

A number of other important points regarding the simulations should also 
be noted. Firstly, in this paper we only examine the 2D case (/3 = 0) where 
rays are confined to the x-z plane. Furthermore, by ensuring that the rays 
remain on the fast- wave branch at all times, we avoid any mode-conversion effects 
as rays pass through the a = c s layer (where fast/slow conversion occurs, see 
Figure 4) . Of course, as numerous works exploring MHD mode conversion in local 
helioseismology have shown (e.g. Spruit and Bogdan, 1992; Cally and Bogdan, 
1993; Cally, Bogdan, and Zweibel, 1994; Bogdan and Cally, 1997; Cally and 
Bogdan, 1997; Cally, 2000, 2007; Crouch and Cally, 2003, 2005; Schunkcr and 
Cally, 2006), mode transmission and conversion between fast and slow magneto- 
acoustic waves indeed occurs as rays of helioscismic interest pass through the 
a = c s cquipartition level and have distinct effects on helioseismic waves that 
should not be ignored. But in our current analysis (and as with actual time- 
distance inversions) we do not directly account for these effects. As a result 
the complexities of the ray-path calculations are greatly reduced. We also note 
that we ignore any finite-wavelength effects and filtering of observations in our 
simulations. 

The computational ray propagation grid extends across the 16 Mm radius 
of the sunspot model in regular 1 Mm spatial increments in the horizontal re- 
direction and down to a depth of 25 Mm in the vertical z-direction, employing a 
much finer grid spacing in the top 1.5 Mm, followed by 1 Mm increments down 
to a depth of 25 Mm. The cutoff height (depth) for all rays propagated in the 
grid was fixed at z = —0.1 Mm, regardless of frequency This computational 
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grid, though not exhaustive, allows us to obtain the desired range of skip dis- 
tances required to replicate the "centre-to-annulus" skip distance geometry (i.e. 
averaging rays from a central point/pixel to a surrounding annulus of different 
sizes to probe varying depths beneath the solar surface) often employed in time- 
distance helioseismology for the derivation of mean travel-time perturbation 
maps (see Gizon and Birch (2005) for a more comprehensive description of this 
process). The 11 standard skip distance bin/annuli (A) sizes usually used for 
these calculations are detailed in Table 1. 

Table 1. The annuli 
(or skip-distances) 
geometries used to bin 
the ray travel-time 
measurements. 



A 


Pupil Size (Mm) 


1 


3.7 - 8.7 


2 


6.2 - 11.2 


3 


8.7 - 14.5 


4 


14.5 - 19.4 


5 


19.4 - 29.3 


6 


26.0 - 35.1 


7 


31.8 - 41.7 


8 


38.4 - 47.5 


9 


44.2 - 54.1 


10 


50.8 - 59.9 


11 


56.6 - 66.7 



4. Results 

4.1. Travel-Time and Skip-Distance Perturbations 

The ray propagation grids were computed for three frequencies, ui = 3.5, 4, and 
5 mHz. Both the phase (i p , associated with the phase velocity) and group (t g , 
associated with the envelope peak of a wave packet as it travels at the group 
velocity) ray travel times were calculated along each ray path for every radial 
grid position (r spot , which is essentially the radial position of the lower turning 
point of the ray) along the sunspot model. In time-distance helioseismology, 
centre-to-annulus travel times are extracted from Gaussian wavelet fits - usually 
represented by a function of the form 

W±(t) = Ac-t 2 ^^ 2 cos[wo(iT* P )], (20) 

(where all parametres are free) - to both the positive and negative time parts 
of the observed cross-correlations (Gizon and Birch, 2005). However, t p is more 
often used in time-distance literature, primarily as a result of difficulties (mainly 
observational noise) associated with fitting to the envelope peak. Furthermore, 
because t p is much more independent of the shape of the wave packet than t g 
(as the shape of the wavepacket depends on (unmodcllcd) mode conversion), 
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we shall also limit our analysis to t p calculations in this paper. We identify the 
phase travel time as 

t P = ^, (21) 

which is consistent with the form of t p described by the Gaussian wavelet. These 
travel times are then subtracted from similar ray travel times calculated using 
the quiet-Sun atmosphere to produce travel-time perturbation (6t p ) profiles. In 
general, travel-time differences are sensitive to sub-surface flows, while mean 
travel times are sensitive to wave-speed perturbations. However, as our model 
does not contain flows, we do not need to distinguish directions along ray paths. 




Figure 6. Travel-time perturbations (6t p ) as a function of skip distance (x) for r spo t = 4, 8, 12, 
and 16 Mm on the sunspot (where r spo t is the radial position of the lower turning point of the 
ray), as calculated for three frequencies: ui = 3.5 (green), u) = 4 (red) and ui = 5 mHz (blue). 

In Figure 6 we see some sample 6r p profiles for r spot = 4, 8, 12, and 16 Mm are 
shown as a function of ray skip distance (x) for ui = 3.5 (green), 4 (red), and 5 
mHz (blue). By and large, there are significant perturbations as we approach the 
centre of the sunspot (i.e. regions associated with stronger surface magnetic field 
strength). The sign of the perturbations appears to remain exclusively negative, 
regardless of position on the sunspot. This means that all rays propagated within 
the simulated sunspot atmosphere are significantly sped up when compared to 
their Model S counterparts. 

Furthermore, in Figure 7 we can see that there are also significant skip- 
distance perturbations (6x) associated with rays that are propagated through the 
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sunspot atmosphere. These calculations are for similar positions and frequencies 
as in Figure 6. The exclusively positive values of 6a; that we can see along the 
sunspot radius indicates that at the same time that these rays are being sped 
up, they are also undertaking a longer journey than their Model S counterparts 
in the process, and as with 6t p , the magnitude of the calculated 6a; appears to 
be closely related to surface magnetic field strength. For both 6r p and 6a; we also 
observe a particular pattern of perturbation associated with each position along 
the sunspot. Whereas the perturbations appear to mainly decrease when we are 
close to spot centre (e.g. r spot = 4, 8 Mm), they appear to increase when further 
away (e.g. r spot = 12, 16 Mm) from spot centre. This is clearly a bi-product of 
both varying field strength and inclination angle of field lines (see Figure 1) as 
we move across the sunspot. Field strength tends to decrease, while field lines 
become more significantly inclined as we move away from centre of the sunspot. 



r spot = 4 Mm r spot = 8 Mm 




10 15 20 25 30 35 10 15 20 25 30 35 

t p (minutes) t p (minutes) 



Figure 7. Skip distance perturbations (5 X ) as a function of phase travel time (t p ) for 
''spot = 4,8, 12, and 16 Mm on the sunspot, calculated for three frequencies^ = 3.5 (green), 
u) = 4 (red), and oj = 5 mHz (blue). 

Also clearly obvious from both Figures 6 and 7 is the presence of a significant 
frequency dependence of both 6t p and 6a; measurements in the sunspot, with the 
magnitudes of the perturbations increasing as the frequency is increased from 
3.5 to 5 mHz. This is particularly evident for rays with short skip distances (i.e. 
surface skimmers with very shallow lower turning points). Frequency dependence 
of travel-time perturbations in active regions has also been observed by both 
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hclioseismic holography (Braun and Birch, 2006) and time-distance hclioseis- 
mology (Couvidat and Rajaguru, 2007). We shall discuss the importance of 
these observations in greater detail in the upcoming sections. Cally (2007) also 
observed a similar behaviour when modelling rays in inclined fields and described 
several related but distinct effects that strong magnetic fields appear to have on 
seismic waves, with an important "dual effect" that the magnetic field has on 
individual ray paths (that is, increasing their skip distances while at the same 
time, speeding them up considerably) being one of these effects. 

A comparison between rays propagated inside the sunspot model with rays 
propagated in the quiet-Sun clearly reveals these effects to the naked eye. All rays 
shown in Figure 8 are initialized at a depth of z^ ot — —2 Mm, with the rays inside 
the sunspot model (solid rays, colours identify frequencies) also being initialized 
at varying positions along the sunspot (r spot = 0, 4, 8, 12, and 16 Mm). While the 
rays propagated inside the Model S atmosphere (dashed rays) are symmetrical 
about their turning points (as expected), strong asymmetries (at both turning 
points) are associated with the same rays when initiated inside the sunspot. 
We can clearly see that the rays inside the sunspot (at all three frequencies) 
appear to have undergone a longer skip distance, in a slightly shorter amount 
of time (dots along ray paths indicate one-minute t g intervals), confirming the 
perturbation profiles of Figures 6 and 7. Of course Figure 8 shows a very small 
sample of rays initialized at a given depth, but even so, they are quite clearly 
indicative of the large-scale effects of the magnetic field on ray propagation - 
effects which are more pronounced as we approach the spot centre and in regions 
of significantly inclined magnetic fields. 

4.2. Binned Travel-Time Perturbation Profiles 

The mean ray travel-time perturbations (6t™) for each frequency and grid po- 
sition were calculated and binned into 11 annuli (Ai — An) of various sizes 
(outlined in Table 1). The 6t™ profiles of the bins are shown in Figure 9. Once 
again, we can see the clear frequency dependence of travel-time perturbations 
evident in all bins, with perturbations increasing with increasing frequency as 
before. Also, all 6t™ bins contain negative perturbations as we saw before in 
Figure 6. We also observe that the magnitude of 6r™ decreases as we move away 
from the centre of the sunspot (i.e. decreasing field strength) for the smaller bins 
{e.g. Ai - A 3 ). 

These smaller bins are representative of shallow rays that spend a considerable 
proportion of their journey inside the magnetic field, consistent with the larger 
magnitude of the perturbations seen in these bins. Larger bins (e.g. A4 — An) 
sample rays with much deeper lower turning points, hence a considerable amount 
of the journey undertaken by these rays would be spent in the quiet-Sun Model 
S atmosphere. Therefore the magnitude of the perturbations tends to be smaller 
than that for the smaller bins. However, they are found to increase in magnitude 
as we move away from the centre of the sunspot as rays sample larger areas of 
the magnetic field throughout their journey across the sunspot radius. 

It should be noted that for the smaller bins (particularly for Ai — A 3 ), it 
becomes quite difficult to obtain a sufficient sampling of rays to average near 
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Figure 8. Individual rays propagated through the simulated sunspot (solid rays) and Model 
S (dashed rays) atmospheres, calculated for three frequencies: u> = 3.5 (green), u> = 4 (red), 
and lu = 5 mHz (blue). The top of each frame indicates the initial depth (zboti Mm) and radial 
grid position of the lower turning point of the ray (r sp ot> Mm). 

the centre of the flux tube, even with a very fine grid spacing of Az — —0.025 Mm 
in the very sensitive top 1.5 Mm of the computational grid. As such, we get a 
certain level of rigidity in the 6t™ profiles of these bins. No such restriction is 
encountered when using a pure Model S atmosphere, which tends to suggest 
that strong near-surface magnetic fields are severely restricting the propagation 
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of helioseismic rays with short skip distances (or very shallow lower turning 
points) . 







Figure 9. Binned (mean) travel-time perturbation (St™, minutes) profiles as a function of 
position (r spot , Mm) on the sunspot, calculated for three frequencies: ui = 3.5 (green), u> = 4 
(red), and u> = 5 mHz (blue). Annuli number and sizes are indicated on the top of the frame 
of each bin. 

Although our sunspot model has many of the qualitative features we might 
expect in a real spot, it is nonetheless rather ad hoc, and consequently our time- 
distance results do not warrant detailed comparison with solar observations. 
Nevertheless, it is of interest to qualitatively compare the 6r™ results obtained 
from our simulations to those reported for AR 8243 (18 June 1998) by Couvidat, 
Birch, and Kosovichev (2006). S. Couvidat kindly provided us with the actual 
set of travel time maps used in their analysis. 

To compare the 6r™ profiles as closely as possible, we first compute the 
azimuthal average of the four 6r™ maps presented in Figure 3 of Couvidat, 
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Birch, and Kosovichev (2006) (corresponding to Ai, A3, A 6 and A 9 , noting 
that the travel-times were obtained without a frequency bandpass filter), to 
obtain 6r™ profiles of AR 8243, akin to our artificial 6t™ profiles contained in 
Figure 9. We observe peak (positive) travel-time perturbations of w 0.29 and 
w 0.16 minutes respectively for A x and A 3 in the sunspot umbra, while the 
sign of 6t™ in the sunspot changes for the larger bins, A§ and Ag, with 6r™ 
ranging from w —0.38 to —0.31 minutes respectively. The perturbations for 
all four bins also appear to decrease in the penumbra relative to the umbra. 
In comparison, if we assume a central frequency of 3.5 mHz, the artificial 6r™ 
profiles for the bins produced by our simulations (Figure 9, 3.5 mHz profiles 
indicated by solid green lines) show opposite-in-sign and larger-in-magnitude 
5t™ for both Ai (w —0.7 minutes) and A 3 (w —0.82 minutes), while similar- 
in-sign yet smaller-in-magnitude 6r™ profiles were observed for A 6 (w —0.22 
minutes) and Ag (m —0.05 minutes). When we consider higher frequencies, the 
magnitude of the artificial 6r™ increases with frequency for all four bins, with 
all perturbations being negative in sign. However, the general pattern of the 
artificial 6r™ profiles for all frequencies appears to be similar to the observations 
of Couvidat, Birch, and Kosovichev (2006), with perturbations decreasing with 
increasing radius from the centre of the sunspot. 

While the differences in the magnitudes of 5r™ between our simulations and 
those of Couvidat, Birch, and Kosovichev (2006) (at a given fixed central fre- 
quency) can be explained, to some extent, by magnetic and thermal differences 
between our model and their sunspot, the frequency dependence of 6r™ and 
the sign change of the smaller bins in particular (i.e. positive 6t™ resulting from 
actual time-distance observations, negative 6r™ from the simulations) can not be 
dismissed as easily. Traditionally, positive 6r™ obtained for short skip distances 
in sunspots have been interpreted as representing a region of slower wave-speed 
propagation in the shallow sub-surface layers of the sunspot. However, as we 
briefly noted in the previous section, Braun and Birch (2006) (using helioseismic 
holography) found that, at a given fixed phase speed, travel-time perturba- 
tions within active regions exhibit a strong frequency dependence. Couvidat and 
Rajaguru (2007) confirmed these results using time-distance hclioscismology, 
applying additional frequency bandpass filters (centred at 3, 4 and 4.5 mHz) to 
the standard phase-speed filters used in Couvidat, Birch, and Kosovichev (2006) 
in order to determine the cause of the dark rings of negative 6r™ they detected in 
the travel-time maps (mainly associated with the A 2 and A 3 skip-distance bins) 
of a majority of the sunspots they studied. These rings, which are sensitive to the 
frequency filtering applied, are found to produce significant ring-like structures 
in the inversion results, mimicking regions of increased sound speed. The authors 
conclude that the rings are most likely to be artifacts caused by surface effects, 
probably of magnetic origin. 

In addition to these results, the very recent work undertaken by Braun and 
Birch (2008) (using ridge filters, in addition to the standard phase-speed filters) 
provide strong evidence that the positive perturbations observed arise from the 
pi ridge or beneath it. These positive travel-time shifts were not seen in the 
higher order p-mode data. These results, when considered in conjunction with 
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our artificial 6r™ profiles (and the results contained in in the next section) , pro- 
vide further concrete evidence that positive travel-time perturbations obtained 
for short skip distances are likely to be artifacts or bi-products of the data 
reduction or analysis method used, rather than some actual physical sub-surface 
anomaly below the sunspot. 

4.3. Isolating the Thermal Component of Travel Time Perturbations 

One of the keys to understanding the role played by near-surface magnetic fields 
in local hclioscismology is to be able to isolate it from effects thought to be 
produced by thermal or flow perturbations. The simplest way to isolate such 
effects is to essentially "switch off" the magnetic field when calculating the ray 
paths in the simulations - that is, set a = in the simulated sunspot atmosphere, 
but maintain the modified sound-speed profile obtained (seen in Figure 4). 

The external atmosphere, ray-path simulations and computational grid re- 
main identical to those described previously. The only difference is the resulting 
thermal travel-time perturbations (Sr™') which would then be purely a result of 
what can be referred to as "thermal variations" along the ray path. One can then 
compare the resulting perturbation profiles to those obtained when the magnetic 
field is included in the simulations (i.e. Figure 9) to better understand the role 
of the thermal contributions to the observed 8t™ profiles. Figure 10 shows the 
resulting bins of the thermal component of 6r™* . 

In general, the resulting 6t™* profiles are relatively smooth and all bins clearly 
show exclusively positive travel-time perturbations (compared to exclusively neg- 
ative travel-time perturbations observed in Figure 9), this implies that rays are 
travelling considerably slower than in the Model S atmosphere - a clear contrast 
with simulations where the magnetic field is present. The magnitude of 6r™* is 
also decreasing with increasing radius for the smaller bins (Ai — A4) and vice 
versa for the larger bins (A 5 — An), a similar behaviour to what is observed in 
Figure 9. However, when considering the magnitude of the perturbations between 
Figures 9 and 10, it is clear that thermal perturbations appear to be much smaller 
for a majority of the bins - in fact up to 400% smaller for some frequencies when 
comparing the perturbations in the near-surface regions Ai— A 3 . The magnitude 
of the perturbations become much more comparable when looking at the larger 
bins ( A7 onwards) , and from A 8 onwards 6r™* becomes ever slightly larger than 
the ones we see in Figure 9 for the same bins. Frequency dependence of 6r™* is 
also evident, but only clearly discernible for the first six bins (Ai — A 6 ). 

5. Summary and Discussion 

Whether it be through direct observations, forward modelling, or inversions, 
in order to be able to confidently interpret helioseismic observations and infer- 
ences made in regions of strong magnetic field, the actual physical effects of 
near-surface magnetic fields on ray propagation must be better understood and 
taken into account when analyzing or modelling active region sub-photospheres. 
Our approach here is akin to forward modelling of rays, but in a simulated 
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Figure 10. Binned (mean) thermal travel-time perturbation (St™*, minutes) profiles as a 
function of position (r spo t, Mm) on the sunspot, calculated for three frequencies: ui = 3.5 
(green), w = 4 (red), and ui = 5 mHz (blue). Annuli number and sizes are indicated on the top 
of the frame of each bin. 



sunspot atmosphere based on IVM surface magnetic-field profiles with a peak 
field strength of 3 kG and an external field-free Model S atmosphere used as the 
background or unperturbed medium. The main aim of these simulations was to 
isolate and understand the effects of the wave-speed inhomogeneities produced 
by the magnetic field from those thought to be produced from thermal or flow 
perturbations. 

The magneto-acoustic rays were propagated across the sunspot radius for a 
range of depths to produce a skip distance geometry similar to centre-to-annulus 
cross-covariances used in time-distance helioseismology. The perturbations from 
the Model S atmosphere were calculated for each radial grid position and range of 
frequencies (3.5 — 5 mHz), then binned into 11 different skip-distance geometries 
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of increasing size. A separate, yet similar, set of simulations was then produced 
to isolate the role played by thermal variations inside the sunspot atmosphere 
on the ray skip-distance and travel-time perturbation profiles. This was achieved 
by having the magnetic field switched off in the sunspot model, thus essentially 
maintaining a modified sound-speed structure, but with no calculations of the 
Alfvcn speed. 

These artificial skip-distance and travel-time perturbation profiles, which di- 
rectly account for the effects near-surface magnetic fields and thermal variations 
separately, have provided us with a number of very distinct and interesting 
observations: 

1 . The sunspot magnetic field has a clear and distinct "dual effect" on helioseis- 
mic rays - increasing their skip distances, while at the same time, shortening 
their travel time (compared to similar rays in a Model S atmosphere). Higher 
frequency rays propagated within the magnetic field also tend to undergo a 
more substantial speed up than their non-magnetic counterparts. 

2. There is a clear and significant frequency dependence of both ray skip-distance 
and travel-time perturbations across the simulated sunspot atmosphere. This 
frequency dependence of perturbations was prevalent for all skip-distance 
bins, but particularly so for shallow rays, which sample the near-surface layers 
of the sunspot. 

3. The negative sign of travel-time shifts, along with the general pattern and 
magnitude of these perturbations (i.e. tending to increase with increasing 
magnetic-field strength and inclination) points to more evidence of the sig- 
nificant role played by the sunspot magnetic field. Rays with shorter skip 
distances were seen to experience greater perturbations as a result of spending 
a considerable proportion of their journey within the confines of the magnetic 
field. 

4. With the magnetic field switched off, the simulated travel-time perturbation 
profiles changed sign for all bins (i.e. only positive perturbations were ob- 
served across the sunspot radius, meaning that rays in the thermal model are 
actually slower than their Model S counterparts), and the magnitude of these 
perturbations appeared to be significantly smaller in magnitude (300-400% 
at times) than when the magnetic field is included in the model. This was 
particularly evident for the bins that sample rays in the near-surface layers, 
whereas bins of larger skip distances produce slightly larger perturbations 
than the magnetic model. Frequency dependence of travel-time perturbations 
were also observed, but only for half of the bins. A majority of bins sampling 
larger skip distances did not exhibit this behaviour. 

These observation as a whole tend to suggest that active-region magnetic fields 
play a direct and significant role in sunspot seismology, and it is the interaction of 
the near-surface magnetic field with solar oscillations, rather than purely thermal 
(or sound-speed) perturbations, that is the major cause of observed travel- 
time perturbations in sunspots. (We note here that we are only commenting 
on the interpretation of time-distance results in terms of thermal/sound-speed 
perturbations, and not, for example, in terms of wave-speed perturbations). 
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The frequency dependence of these perturbations is one of the strongest 
indications that the magnetic field is a significant contributor to the travel- 
time shifts. When isolating the thermal component of St p we did observe some 
frequency dependence in a limited number of bins/skip distance geometries, 
certainly not to the extent that we saw when the magnetic field was included. 
Of course in the absence of any perturbations, rays propagated at different 
frequencies will naturally have slightly different upper turning points, this could 
certainly explain a part of a frequency dependence, but this effect combined 
with the (negative) sign and magnitude of the simulated 5r p profiles, along with 
the relatively small (positive) thermal component extracted from the perturba- 
tions, makes it very difficult for one to argue that what we arc seeing in these 
travel-time perturbation profiles is a result of a sub-surface flow or sound-speed 
perturbation, as has been traditionally interpreted in time-distance literature. 

Instead, these observations indicate that strong near-surface magnetic fields 
may be seriously altering the magnitude and lateral extent of sound-speed in- 
versions made by time-distance helioseismology. This is because standard time- 
distance observations (e.g. Couvidat, Birch, and Kosovichev (2006), see Section 
4.2) show St™ maps derived from the averaged cross-correlations shifting from 
positive values for the first couple of bins (usually Ai — A3), to negative ones for 
the remainder of the bins. Traditionally, positive perturbations result in regions 
of decreased sound speed in inversions, while negative perturbations result in 
regions of enhanced sound speed. But we have clearly seen from our forward 
modelling that the inclusion of the magnetic field in the near surface layers 
consistently results in negative values for all bins of St™. This implies that any 
inversion of time-distance data that does not account for surface magnetic field 
effects will be significantly contaminated in the shallower layers of the sunspot 
(i.e. down to a depth of a few Mm below the surface), in strong agreement with 
the conclusions of Couvidat and Rajaguru (2007). Hence it is almost certain 
from these simulations that the two-structure sunspot sound speed profile, i.e. 
region of decreased sound speed immediately below the sunspot (corresponding 
to positive St™), is most likely an artifact due to surface effects, instead of 
thermal perturbations. Deeper sound speed profiles do not appear to be affected 
as much, given the sign and magnitude of the simulated St™ for the larger bins 
are comparable to actual time-distance calculations, as expected, given the flux 
tube becomes gas-pressure dominated at such depths. 

Of course, we must bear in mind that some of our assumptions outlined earlier 
(e.g. 2D treatment of rays, our choice of ray cutoff height in the atmosphere, the 
fact that we are not directly accounting for mode conversion, even the form of the 
surface magnetic field and background model in general etc.), can certainly alter 
our results quantitatively in one manner or another. Indeed it would certainly 
be interesting and worthwhile to conduct a full 3D simulation (i.e vary the 
shooting angle (3 around the sunspot) and also test the ray propagation code 
with other sunspot and quiet-Sun models in the future. But in any case, it would 
be surprising, given the self-consistency of our current results, if our qualitative 
conclusions were changed as a result. 
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